{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import\n",
    "from collections import defaultdict\n",
    "from functools import lru_cache\n",
    "from tqdm import tqdm\n",
    "import qiskit.circuit.random\n",
    "import scipy.io as sio\n",
    "import random\n",
    "from qiskit.ignis.mitigation.measurement import complete_meas_cal, CompleteMeasFitter\n",
    "import numpy as np\n",
    "import scipy.linalg as la\n",
    "from qiskit.quantum_info import hellinger_fidelity\n",
    "from qiskit.result import LocalReadoutMitigator\n",
    "from qiskit.visualization import plot_histogram\n",
    "from qiskit import QuantumCircuit\n",
    "import matplotlib.pyplot as plt\n",
    "import jax\n",
    "from qiskit.result import Counts\n",
    "from jax import numpy as jnp\n",
    "from jax import vmap\n",
    "from qiskit_aer import AerSimulator\n",
    "# from jax import random\n",
    "from qiskit_aer.noise import NoiseModel, ReadoutError\n",
    "import random\n",
    "import math\n",
    "import time\n",
    "import psutil\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # def cartesian_coord(arr1, arr2):\n",
    "# #     grid = np.meshgrid(arr2, arr1)\n",
    "# #     coord_list = [entry.ravel() for entry in grid]\n",
    "# #     points = np.vstack(coord_list).T\n",
    "# #     return points\n",
    "\n",
    "# # @jax.jit\n",
    "# def kron_basis(arr1, arr2):\n",
    "#     grid = jnp.meshgrid(arr2, arr1)\n",
    "#     return grid[1].ravel() << 1 | grid[0].ravel()\n",
    "\n",
    "\n",
    "# basis_1q = jnp.array([[0, 1], [1, 0]])\n",
    "    \n",
    "# @jax.jit\n",
    "# def process_sample(local_basis, qubit, now_basis, now_values, all_local_vecs, threshold): \n",
    "#     next_basis = kron_basis(now_basis, basis_1q[local_basis])\n",
    "#     next_values = jnp.kron(now_values, all_local_vecs[qubit][local_basis])\n",
    "#     filter = jnp.logical_or(next_values > threshold, next_values < -threshold)\n",
    "#     return next_basis, next_values, filter\n",
    "\n",
    "# def sw_rm(qnum, stats_counts, meas_mats_inv, threshold=1e-10):\n",
    "#     all_local_vecs = np.zeros(shape=(qnum, 2, 2))\n",
    "\n",
    "#     for qubit in range(qnum):\n",
    "#         for local_basis in (0, 1):\n",
    "#             local_vec = np.zeros(2)\n",
    "#             local_vec[local_basis] = 1\n",
    "#             local_vec = meas_mats_inv[qubit] @ local_vec\n",
    "#             all_local_vecs[qubit][local_basis] = local_vec\n",
    "\n",
    "#     all_local_vecs = jnp.array(all_local_vecs)\n",
    "\n",
    "\n",
    "#     rm_prob = defaultdict(float)\n",
    "#     # for basis, count in tqdm(stats_counts.items()):\n",
    "#     for basis, count in stats_counts.items():\n",
    "#         basis = jnp.array([int(c) for c in basis])\n",
    "\n",
    "#         now_basis = basis_1q[basis[0]]\n",
    "#         now_values = all_local_vecs[0][basis[0]] * count\n",
    "        \n",
    "#         for qubit in range(1, qnum):\n",
    "#             next_basis, next_values, filter = process_sample(basis[qubit], qubit, now_basis, now_values, all_local_vecs, threshold)\n",
    "#             now_basis, now_values = next_basis[filter], next_values[filter]\n",
    "\n",
    "#         for basis, value in zip(now_basis, now_values):\n",
    "#             basis = int(basis)\n",
    "#             # print(basis, value)\n",
    "#             rm_prob[basis] += value\n",
    "\n",
    "#     sum_prob = sum(rm_prob.values())\n",
    "#     rm_prob = {\n",
    "#         basis: value / sum_prob\n",
    "#         for basis, value in rm_prob.items()\n",
    "#     }\n",
    "#     return rm_prob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "# def cartesian_coord(arr1, arr2):\n",
    "#     grid = np.meshgrid(arr2, arr1)\n",
    "#     coord_list = [entry.ravel() for entry in grid]\n",
    "#     points = np.vstack(coord_list).T\n",
    "#     return points\n",
    "\n",
    "\n",
    "def kron_basis(arr1, arr2):\n",
    "    grid = np.meshgrid(arr2, arr1)\n",
    "    return grid[1].ravel() << 1 | grid[0].ravel()\n",
    "\n",
    "def sw_rm(qnum, stats_counts, meas_mats_inv, threshold=1e-10):\n",
    "    all_local_vecs = np.zeros(shape=(qnum, 2, 2))\n",
    "\n",
    "    for qubit in range(qnum):\n",
    "        for local_basis in (0, 1):\n",
    "            local_vec = np.zeros(2)\n",
    "            local_vec[local_basis] = 1\n",
    "            local_vec = meas_mats_inv[qubit] @ local_vec\n",
    "            all_local_vecs[qubit][local_basis] = local_vec\n",
    "\n",
    "    basis_1q = jnp.array([[0, 1], [1, 0]])\n",
    "\n",
    "    rm_prob = defaultdict(float)\n",
    "    # for basis, count in tqdm(stats_counts.items()):\n",
    "    for basis, count in stats_counts.items():\n",
    "        basis = [int(c) for c in basis]\n",
    "\n",
    "        now_basis = basis_1q[basis[0]]\n",
    "        now_values = all_local_vecs[0][basis[0]] * count\n",
    "\n",
    "        for qubit in range(1, qnum):\n",
    "            next_basis = kron_basis(now_basis, basis_1q[basis[0]])\n",
    "            next_values = np.kron(now_values, all_local_vecs[qubit][basis[qubit]])\n",
    "            \n",
    "            filter = np.logical_or(next_values > threshold, next_values < -threshold)\n",
    "            now_basis = next_basis[filter]\n",
    "            now_values = next_values[filter]\n",
    "\n",
    "        for basis, value in zip(now_basis, now_values):\n",
    "            rm_prob[basis] += value\n",
    "\n",
    "    sum_prob = sum(rm_prob.values())\n",
    "    rm_prob = {\n",
    "        basis: value / sum_prob\n",
    "        for basis, value in rm_prob.items()\n",
    "    }\n",
    "    return rm_prob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\n",
      "10\n",
      "11\n",
      "12\n",
      "13\n",
      "14\n",
      "15\n",
      "16\n",
      "17\n",
      "18\n",
      "19\n",
      "20\n",
      "21\n",
      "22\n",
      "23\n",
      "24\n",
      "25\n",
      "26\n",
      "27\n",
      "28\n",
      "29\n",
      "30\n",
      "31\n",
      "32\n",
      "33\n",
      "34\n",
      "35\n",
      "36\n",
      "37\n",
      "38\n",
      "39\n",
      "40\n",
      "41\n",
      "42\n",
      "43\n",
      "44\n",
      "45\n",
      "46\n",
      "47\n",
      "48\n",
      "49\n",
      "50\n",
      "51\n",
      "52\n",
      "53\n",
      "54\n",
      "55\n",
      "56\n",
      "57\n",
      "58\n",
      "59\n",
      "60\n",
      "61\n",
      "62\n",
      "63\n",
      "64\n",
      "65\n",
      "66\n",
      "67\n",
      "68\n",
      "69\n"
     ]
    }
   ],
   "source": [
    "def ghz(n_qubits):\n",
    "    cir = QuantumCircuit(n_qubits)\n",
    "    cir.h(0)\n",
    "    for i in range(n_qubits - 1):\n",
    "        cir.cx(i, i + 1)\n",
    "    cir.measure_all()\n",
    "    return cir\n",
    "\n",
    "time_qiskit = []\n",
    "time_sw = []\n",
    "space_qiskit = []\n",
    "space_sw = []\n",
    "\n",
    "for n_qubits in range(5, 70):\n",
    "    print(n_qubits)\n",
    "    measure_fids = np.random.random(size=(n_qubits, 2))\n",
    "    measure_fids = np.abs(measure_fids) / 10 + .9\n",
    "    noise_model = NoiseModel()\n",
    "\n",
    "    for i in range(n_qubits):\n",
    "        re = ReadoutError([[measure_fids[i][0], 1 - measure_fids[i][0]],\n",
    "                           [1 - measure_fids[i][1], measure_fids[i][1]]])\n",
    "        noise_model.add_readout_error(re, qubits=[i])\n",
    "    simulator = AerSimulator(noise_model=noise_model)\n",
    "    noise_free_simulator = AerSimulator()\n",
    "    \n",
    "    n_samples = 3000\n",
    "    before_rm_counts = simulator.run(\n",
    "        ghz(n_qubits), shots=n_samples).result().get_counts()\n",
    "    \n",
    "    meas_mats, meas_mats_inv = [], []\n",
    "    for qubit in range(n_qubits):\n",
    "        meas_mat = np.array([[\n",
    "            measure_fids[qubit][0], 1-measure_fids[qubit][1]],\n",
    "            [1-measure_fids[qubit][0], measure_fids[qubit][1]]\n",
    "        ])\n",
    "        meas_mats.append(meas_mat)\n",
    "        meas_mats_inv.append(np.linalg.inv(meas_mat))\n",
    "\n",
    "    # if n_qubits < 16:\n",
    "    #     mit = LocalReadoutMitigator(meas_mats, list(range(n_qubits)))\n",
    "        \n",
    "    #     start_time = time.time()\n",
    "    #     qiskit_rm_prob = mit.quasi_probabilities(before_rm_counts)\n",
    "    #     qiskit_rm_prob = {\n",
    "    #         basis: value\n",
    "    #         for basis, value in qiskit_rm_prob.items()\n",
    "    #         if value > 1e-3\n",
    "    #     }\n",
    "    #     time_qiskit.append(time.time() - start_time)\n",
    "        \n",
    "    start_time = time.time()\n",
    "    rm_prob = sw_rm(n_qubits, before_rm_counts, meas_mats_inv, threshold=n_samples * 1e-5)\n",
    "    time_sw.append(time.time() - start_time)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# time_sw = [time_sw[i] for i in range(70)]\n",
    "plt.plot(list(range(len(time_qiskit))), time_qiskit)\n",
    "plt.plot(list(range(len(time_sw))), time_sw)\n",
    "plt.xlabel(\"Qubit\")\n",
    "plt.ylabel(\"Time (s)\")\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAiIElEQVR4nO3dd5hU5d3G8e+PpUlbpAjSRAVFVKSsYldswYJirKigRiXmjVFjzPuqKcZETUyMLViCvURAsQEWREUUC1KkN6mCUkSKdLb83j+eIVJ2Z5dlZs/M2ftzXXMtc2Zn5l4d9uac55znMXdHRESkJFWiDiAiIplNRSEiIkmpKEREJCkVhYiIJKWiEBGRpKpGHSAdGjVq5K1bt446hohI1hg/fvwKd29c3GOxLIrWrVszbty4qGOIiGQNM1tY0mM69CQiIknFqijMrIeZ9V+zZk3UUUREYiNWReHuQ929b25ubtRRRERiI1ZFISIiqaeiEBGRpGJVFBqjEBFJvVgVhcYoRERSL1ZFISJSKeRvgi8eh5XzK+TtYnnBnYhIrH32L/jgTrAc6HARHHcTNGqbtreL1R6FxihEJPbWr4DRD8L+J0PXa2Haa9DvcBj8M1g2PS1vGaui0BiFiMTeqL9D/gbo/jfofjfcOAWOuQFmD4fHjoW1S1P+ljr0JCKSLb6fC+OehM59oPEBYVudxnDqHaEs5o+Cuk1T/rax2qMQEYm1D/4COdXhxFt2fqxWAzj43LS8rYpCRCQbLB4fxiOO/lVa9hqSiVVRaDBbRDJWwWZYMBqKikr+nqJCeOF8+FcejH82PAfAHUb8EWo3DkVRwWJVFBrMFpGM4w7TXoeHj4BnzoT3bi/5ez/6B8wZAV4EQ6+HBzrAJw/B1Fdg4Wg44f+gRt0Ki76VBrNFRNJl8TgY/jtY9Dns1R4OOhs+fQgat4NOl27/vfM/gg//Bof1gp6PwrwPYfT9MOIP4fEG+0OXKyr6JwBUFCIiqbfpB3jrZpg8CGrvBT0ehE69w57CC+fB0BugwX6wz1Hh+9d9B69cHS6aO+NeMIP9u4XbNxNg3FOhQHKqRfLjxOrQk4hISk1/Ax47Dia/HA4hlcXymfB4N5gyGI67Ga6fEPYEquSEX/QXPAP1W8GgS2HVwjBm8Vpf2LQmPFajzvav17wznNMPWh+T4h+u7FQUIiI7KiqC9/8CL/UJ8ym9ejW81Dv8yz+Z6W/AEyeHX/qXD4GT/7DzmEKtBnDJICgqgAEXw8g7Ye4HcPo90OTg9P1MuyFWRaGznkSkzIoKYcuGnbdvWgMDe8HH90Kny+Dm2XDKHeHK50e6hlNUi3utEbeHYmncDvqOgtbHlvzejdqGvYfvZsHH/4RDzoPOl6fsR0s187LuTmWRvLw8HzduXNQxRCRTrVoIz50Dq7+GpodAy67hVq8ZDLkeVs0PU2QcfnUYLwBYPgNe/wV8+2X43pzqsHktbFkPG1fBhhXQ5cqwZ1C1RtlyfPkCTH01lEbNemn7ccvCzMa7e16xj6koRKRSWbUAnjkLNv8Qxg6+nRjOTspfHx6v1QgufLb4PYLCAvj0QZgxFKrVgup1oHrtcNvvRDj0/Ir7OVIsWVHorCcRqTxWzoNnesCWddBnCDTrGLYXFsDyabB0aviFn9u8+OfnVIXjfhNulYiKQkQqh+/nhj2Jgk1w+VDYu8OPj+VUhb0PCzfZSawGs0VEirVkcrgqunDzziUhpdIehYhkpqJCWLMYqlQN1x/kVAsDyNVq/TjAXJpVC2Hk3YkL3xrD5cOgSfv05o4hFYWIZJ5VC+HlK+DbCTs/tvdhcNYD4UK0kqxfAR/dG9ZusCpwzPVw7K9hjz3TlTjWYlUUZtYD6NGmTZuoo4hIec16B177eZju4id3hzOLivKhMD8MQn/xRLioreu10O1321/JvGQSTHgOJg0Mq8B1ugxOuKXkwWkpE50eKyKZobAARt4Fo++DpofChc+F+ZB2tGkNvHdH2FvIbRnKZP13MOHZUBQ5NaD9OXD8zdD4wIr/ObKUTo8Vkczy/VxYsyj80t+4Onyd/Q4s/CRcoXz6PVBtj+KfWzMXzroPOlwUJtd7qXfY3uQQOP0f0OECHWJKMRWFiFScdcvDAjyTBuz8WI160PMx6NirbK/Vqiv8/COYMSTseTTrVPZBbtklKgoRKZ+CLbBkIjTPgyqlnGlfWBAOFX1wVxg7OPbX0ObUsHewR/3wtXqdXf9FX7V6Vl8NnS1UFCKy67ash0G9Ye770OJwOP3vxZ+F5B4W5Hn3d7B0Cux/Ujg81EgnnGQTFYWI7JqNq+DFi2Dx2HDm0dRX4fGTwhlGJ98OdRrDD0vC4aUvX4CVc6FeizA4fdDZOjyUhVQUIlJ265bD8+eG6bHPfxoO7hlOUR11D4x5DKYPgRZdYN4o8ELY55hw9lH7nlC9VtTppZxUFCJSNqsWwvM9Ye1SuPSlcBgJwvTYP7krnK00/DZYMStc4NapNzTcP9LIkhoqChEp2cbVYYxh3siwt1CUD33egJZH7Py9jQ+AywZXeERJPxWFSFy4p+b4f8EW+KJ/WMnt2wnhCunqdWDf4+GkP2iupEpIRSESB5vXwdPdw6mqPR4o/+vM/wje/A2smA3Nu8BxN4dDTC3ywqR8UinFqig015NUWm//Xzj9dOmUsDLbrl5bsG45vPv7MMtq/X3g0sHQ9tT0ZJWso7meRLLdtNfCTKvH3BimwFgxG37xKeS2SP68okL4ZnyYOmPsE7BlAxx7Y1i9raTpMyS2NNeTSFytWRzmO2reBU76fZg/6dFj4bVrw1KfO14xXbA5rPc8ezjMeQ82rgTLCYeXfnJ3GJAW2YGKQiRbFRXCqz8PX3/6eBhDaLBfmFBvyHXw+cNw9K9+/P6Fn4VSWTELajWEtqeFw0ttTtYkepKUikIkW33yICwcDec8sv31Cp0uC4eT3v8z7HdiOAT13p9g/DOQ2wp6DQwlUSUnouCSbVQUItlo8biwdkP7ntDxku0fM4MeD8GjR8FLfcIZURtWwFHXQbfboHrtSCJL9iplykcRySiF+WGJz6fPgDpNw6mwxV07Ubth2NNYOQ/qNYNrRoarp1USUg7aoxDJFt9MgCG/gmVTw57E6X9PPrbQ9hS4cQrUbQY5+qsu5adPj0im27wOPvwrfP4I1GkCF78I7c4s23Prt0pvNqkUVBQimaqwACa+ACPvhnXLoMsVcModYaEfkQqkohCJ0g9LYP13UL8l1Kwfxhvc4asRYcnQ72ZAy65w0QvFT8QnUgFUFCJRWfQFPHs2FGwM96vXDYVRpSosnRyuidBiP5IBVBQiqbJxFSyZDEsmhVv+Bjj1L8Uv+7l8BvznAqi3d7ii+odvYfWicGX1+hVhoLrLlWFNaJGIqShEdteSSfDKNeGK563qtYAt6+DxbnDuY9sPPq9aGFaJq1oTer8Ge7au8Mgiu0JFIbI7vh4T9gxq1IVT/gRNO8Deh0HtRmEP4aXeMPASOPamsOewYWUoifwNcOXbKgnJCioKkfKaNwoG9IK6TcOqb/Vbbv94/ZZw5Tvw9m9h9H1hEaCNq8Nhpj6vQ5ODo0gtsssyvijMbD/gd0Cuu+/iJPsiaTJ7OAxKrAnd+3Wo26T476tWE87+V1hQ6K3fQlEB9BoArY6s0LgiuyOSKTzM7CkzW25mU3fY3t3MZpnZHDO7BcDd57n7VVHkFNmJO0x+ORxO2usguOLNkktiW10uh5+PgivfggN+kv6cIikU1VxPzwDdt91gZjnAw8DpQHugl5lpcV7JHF+PgWd7wKtXQ4vD4fIhUKtB2Z+/10Hak5CsFMmhJ3f/yMxa77D5CGCOu88DMLOBwDnA9LK8ppn1BfoCtGqlaQskhZZOgQ/uDFN3194rcerqFVC1RtTJRCpEJo1RNAcWbXN/MdDVzBoCdwGdzOxWd/9rcU929/5AfwhLoaY7rMScO8wfBV88DjOHQc1cOPl26PpzzcAqlU4mFUWx3P174Nqoc0glsekHmDQwrCG9Yhbs0QCO/21Yy0FzLEkllUlF8Q2w7fmFLRLbyszMegA92rQp5kpYkdLMfhcGXxkulGvWGXo+BgefG85cEqnEMqkoxgJtzWxfQkFcDFyS/Cnbc/ehwNC8vLxr0pBP4mzjqrDOdP1WcHY/aNEl6kQiGSOq02MHAJ8BB5rZYjO7yt0LgOuA4cAM4CV3nxZFPqmERvwxzLHU81GVhMgOojrrqVcJ298C3irv6+rQk5TL/I9gwnNwzA3QrGPUaUQyTqzWzHb3oe7eNzc3N+ooki3yN8LQG2DPfeGEW6JOI5KRMmmMQqTijboHVs6DPkOgeq2o04hkpFjtUYjskiWT4ZOHoNNlsN8JUacRyVixKgoz62Fm/desWRN1FMl0m9fCkF9BrYZhcSERKVGsikJjFFKqwgIY9zQ81BmWTIQz7921+ZpEKiGNUUj8bFgJa5dA7cZhj6FKTpiS46sRMOIP8N1MaHVUmO67RV7UaUUynopC4qGoKMzNNOFZmDEMivLDdqsSyqJ6bVi1ABrsBxe9AO3OArNII4tki1gVha6jqCTcYdNqWLss7Dl8Mw4mPA+rF8Iee8LhV0PLw2H997D+O1i/HDZ8D0f+D3S5EqpWj/onEMkq5h6/iVbz8vJ83LhxUceQVFs5P6xBveIrKNi0/WOtjwtTf7c7S3MziZSDmY1392KPxcZqj0JizB3e/A2sXABHXAN1moa1qus2DRfL5TaPOqFIbKkoJDvMGAJz34fu98CRmnVepCLF6vRYianNa+HtW6DpoWH8QUQqlPYoJPONugfWfgsXPgc5+siKVLRY7VHoyuwYWjYdPnsEOl8ezmQSkQoXq6LQldkx4w5v3hTWqz7lT1GnEam0tB8vmWvSAPj6s7DinKbZEIlMrPYoJEaWz4Dht0GLI6DjpVGnEanUVBSSeb6bDc+eDTk14NzHoIo+piJR0t9AySzfz4VnewAOlw+FhvtHnUik0otVUeispyyxeR0sGA2bfth++8r5oSSK8sOKc40PiCafiGwnVoPZ7j4UGJqXl3dN1FmkBIUFMPCSMNOrVYEmh4Qpv5t3hg/ugvwNYU+iSfuok4pIQqyKQrLAiD+GkjjxVvCicFbTl8/DF/8Op8H2GRKuwBaRjKGikIozcQB8/jB0vRZOvOXH7YX5sHQK1GsWJvkTkYxSalGYWU3gLOA4oBmwEZgKvOnu09IbT2Jj8XgYekOYDvy0O7d/LKdaOPQkIhkpaVGY2R2EkvgQGAMsB2oCBwB/S5TIb9x9cppzSjZbuxQGXQp1m8AFz4ZiEJGsUdoexRfufnsJj91nZnsBrVKcSeJk7TIY1Bs2rYGrRkDthlEnEpFdlLQo3P3NHbeZWRWgjrv/4O7LCXsZIj9yh8VjYcy/Yfob4IVw/lPQ9JCok4lIOZRpMNvMXgSuBQqBsUA9M3vQ3f+RznC7SmtmZ4Bpr8Po+2HJRKhRL6xGd/jVunBOJIuV9YK79u7+A9ATeBvYF+idrlDlpdljIzZ9CLx8OeRvhDP/CTfNgO5/VUmIZLmynh5bzcyqEYqin7vnm5mnL5ZknRVz4PX/geZ5cOXbULV61IlEJEXKukfxb2ABUBv4yMz2AX5I+gypPLash5d6h7OZLnxWJSESM2UqCnd/yN2bu/sZ7u7A10C39EaTjDLrnTAP05z3wmD1Vu4w7NdhWvDzn4TcFtFlFJG0SFoUZnZZ4iyn7XhQYGb7m9mx6YsnGeOzfjD/I3jhPHi+JyyZFLaPewomD4Jut8H+J0UaUUTSo7QxiobAl2Y2HhgPfEe44K4NcAKwAril5KdLLKxfAQs/gWNugLp7w6i/w79PgIN6wOx3oM2pcNzNUacUkTQp7TqKB82sH3AScAzQgTCFxwygt7t/nf6IErmZb4YJ/A45H/buAIf1gk8egM8fhTpN4af9tbiQSIyVetaTuxcCIxI3qYxmDIH6+/w4q+se9eGUP8GRvwxThWs9a5FY0z8DJbmNq2HeKGh/Npht/1idxpqSQ6QSUFFIcrOHhxXnDjon6iQiEpFYFYWWQk2DGUOgbjNo3iXqJCISkTIVhZk1MbMnzeztxP32ZnZVeqPtOk3hkWKb14XrJg46S4PVIpVYWf/2PwMMJyxcBDAbuDENeSSTzHkPCjbBQWdHnUREIlTWomjk7i8BRQDuXkCYSVbibMYQqNUI9jk66iQiEqGyFsV6M2sIOICZHQloICDO8jeFgex2Z0KVnKjTiEiEyjp77E3AEGB/M/sEaAycn7ZUEr15H8KWdeG0WBGp1MpUFO4+wcxOAA4EDJjl7vlpTSbRmjEEauRC6+OjTiIiESvrCnc5wBlA68RzTjMz3P2+NGaTqBTmh2k7DjxdU4aLSJkPPQ0FNgFTSAxoSwxtWAlTXobxz8Km1XBwz6gTiUgGKGtRtHD3DmlNItH5egyMfQKmvwGFm6FZZzi7HxzQPepkIpIByloUb5vZae7+blrTSGq5h1lfk5219Plj8M7/Qc1c6HI5dO7z4+R/IiKUvSg+B15LLGKUTxjQdnevl7ZksvtG/R0+fxhOuws6Xbb9pH7u8MGd8PG90O6sMFV49drRZRWRjFXW6yjuA44Carl7PXevq5LIcGuXwej7QyEMuQ5e+CmsTiwfUlgAQ28IJdG5D1zwrEpCREpU1j2KRcDUxHrZkg0+/meY9fUXn8DcD2DE7fDIUWEdiXkfwsxhYVW6k36/8/ThIiLbKGtRzAM+TEwKuHnrRp0em6FWL4LxT4fDTQ33D7e2p4W9iLcSS5Z2vweOvDbanCKSFcpaFPMTt+qJm2SyUfcABsf/74/b9twHer8GUwZDjbpwoM5oEpGyKeuV2XekO0hJzKw28AiwBfjQ3f8TVZassGIOTHwRjugLuc23f8wMOlwQTS4RyVpJB7PNrF/i61AzG7LjrbxvamZPmdlyM5u6w/buZjbLzOaY2S2JzT8FBrv7NYAmHirNh3dD1Rpw3E1RJxGRmChtj6IPcB1wb4rf9xmgH/Dc1g2JaUIeBk4FFgNjE2XUgnBFOGhq8+SWToWpr8CxN0GdvaJOIyIxUVpRzAVw91GpfFN3/8jMWu+w+QhgjrvPAzCzgcA5hNJoAUwkZku3ptzIu8JEfsdcH3USEYmR0oqisZmVeAwjxWc9NSechrvVYqAr8BDQz8zOJMw5VSwz6wv0BWjVqlUKY2WozWvDHsTyabBsGiybDos+h26/hz32jDqdiMRIaUWRA9QhXIkdCXdfD1xZhu/rD/QHyMvLi/f1HpvWwMNdYe2ScL9GLjRpD0dfD0f9MtpsIhI7pRXFEnf/c4UkgW+Altvcb5HYJjuaMSyURI+HoM3JUK+5LpoTkbQp7Zh/Rf72GQu0NbN9zaw6cDFhVb0yM7MeZtZ/zZqYr9I6dTDs2TpMv5HbQiUhImlVWlGcnI43NbMBwGfAgWa22MyucvcCwhlWw4EZwEvuPm1XXtfdh7p739zc3NSHzhTrlsO8UXDIeSoIEakQSQ89ufvKdLypu/cqYftbwFvpeM/YmPY6eCEcoiXLRaRixOp000px6GnqYNirfRi8FhGpALEqitgfelr9NSwaA4dqb0JEKk6siiL2pr4Svh5yXrQ5RKRSUVFkkymvQIvDwxlPIiIVJFZFEesxiuUzYdkUDWKLSIWLVVHEeoxi6mCwKnDwuVEnEZFKJlZFEVvuYcGh1sdB3SZRpxGRSkZFkQ2+nQCr5utsJxGJRKyKIrZjFFNegSrV4KAeUScRkUooVkURyzGKJZNg3FPQ7gxNHy4ikYhVUcTO+hUw8FKo1QDOSPUigyIiZVPaNOMSlcJ8ePmKMAngz97R0qYiEhkVRaZ69w+w4GPo+Rg07xx1GhGpxGJ16Ck2g9kTX4Qxj0LXX0DHYifaFRGpMLEqilgMZn/7JQy9MVwzcdpfok4jIhKvooiFD++BmvXggmcgp1rUaUREVBQZZf33MGcEHHYx1G4UdRoREUBFkVmmvQpFBdDh4qiTiIj8l4oik0weBHsdDE0PiTqJiMh/xaoosvqsp+/nwuKxcNhFUScREdlOrIoiq896mjwIMDj0gqiTiIhsJ1ZFkbXcQ1HsezzUaxZ1GhGR7agoMsGiL2DVgnC2k4hIhlFRZILJg6DqHppGXEQykooiagVbwmmx7c6EGnWjTiMishMVRdTmjICNq3TYSUQylooiapMGQu3GsF+3qJOIiBQrVkWRdddRrP4aZr8Dh5wPOZrxXUQyU6yKIiuuo1j9NXzaD548DR7oAF4EnS6NOpWISIn0z9iKkr8RXrwQ5n8U7jc5FLrdBu17QuMDIo0mIpKMiqKijH0ylMQJt0CHC6Hh/lEnEhEpExVFRdiyHkbfD/udCN1ujTqNiMguidUYRcb64nHYsAJOvC3qJCIiu0xFkW6b18InD0KbU6BV16jTiIjsMhVFuo15DDauDAPXIiJZSEWRTpvWwKf/ggNOh+Zdok4jIlIuKop0+uyRUBYawBaRLKaiSJcNK+HzR8KMsHsfFnUaEZFyU1Gky2cPh4HsE7U3ISLZLVZFkTFzPeVvgrFPhL2JJgdHm0VEZDfFqigyZq6nGUNh02o4/Opoc4iIpECsiiJjfPkc1N8HWh8XdRIRkd2moki1lfPCnE6de0MV/ecVkeyn32Sp9uULYFWgo6YOF5F4UFGkUmEBTHwR2pwK9ZpFnUZEJCVUFKk05z1YuyQcdhIRiQkVRSp9+XxY//qA7lEnERFJGRVFqqxdBrPehsN6QU61qNOIiKSMiiJVJr0IXgid+0SdREQkpVQUqeAOE56HVkdBo7ZRpxERSSkVRSos/BRWztXehIjEkooiFSYNgOp1of05UScREUk5FcXuKiyAWW/BAT+B6rWjTiMiknIqit21aAxs+B7anRl1EhGRtFBR7K6ZwyCnBrQ9NeokIiJpkfFFYWb7mdmTZjY46iw7cQ9Fsd+JUKNu1GlERNIirUVhZk+Z2XIzm7rD9u5mNsvM5pjZLclew93nuftV6cxZbkunwOqvddhJRGKtappf/xmgH/Dc1g1mlgM8DJwKLAbGmtkQIAf46w7P/5m7L09zxvKbOSzMFHvgGVEnERFJm7QWhbt/ZGatd9h8BDDH3ecBmNlA4Bx3/ytwVnnfy8z6An0BWrVqVd6X2TUz34SWR0KdxhXzfiIiEYhijKI5sGib+4sT24plZg3N7DGgk5ndWtL3uXt/d89z97zGjSvgF/fK+bBsqg47iUjspfvQ025z9++Ba6POsZOZb4avB5V7J0hEJCtEsUfxDdBym/stEtt2m5n1MLP+a9asScXLJTdzGDQ5FPZsnf73EhGJUBRFMRZoa2b7mll14GJgSCpe2N2Hunvf3NzcVLxcydYth68/12EnEakU0n167ADgM+BAM1tsZle5ewFwHTAcmAG85O7T0pkj5Wa9DbgOO4lIpZDus556lbD9LeCtdL53Ws0cBvX3gSaHRJ1ERCTtMv7K7F1RIWMUm9fCvA+h3Vlglr73ERHJELEqigoZo/j4PijcAgf1SN97iIhkkFgVRdp99giMvi8sUNTqyKjTiIhUiFgVRVoPPU0cAMNvhYPOhrMe0GEnEak0YlUUaTv0NPMteOOXYZbY856AKjmpfX0RkQwWq6JIi/kfw8tXQLOOcNF/oGqNqBOJiFQoFUUyc0fCgF7QYF+4dDDUqBN1IhGRCqeiKMnEAfCf86F+S+j9GtRqEHUiEZFIxKooUjKY7Q6j/gGvXwv7HA0/ewfqNUtdSBGRLBOrotjtwezCfBh6PYy8EzpcDJe+AjXTPG+UiEiGy/hpxitMwRYYeAnMGQHH3Qwn/V6nwIqIoKL4UU41aNgmzAibd2XUaUREMoaKYiszOP1vUacQEck4sRqjqNCFi0REKolYFUWFLVwkIlKJxKooREQk9VQUIiKSlIpCRESSUlGIiEhSsSoKnfUkIpJ6sSoKnfUkIpJ65u5RZ0g5M/sOWFjOpzcCVqQwTkXK5uyg/FHK5uyQ3fkzJfs+7t64uAdiWRS7w8zGuXte1DnKI5uzg/JHKZuzQ3bnz4bssTr0JCIiqaeiEBGRpFQUO+sfdYDdkM3ZQfmjlM3ZIbvzZ3x2jVGIiEhS2qMQEZGkVBQiIpKUiiLBzLqb2Swzm2Nmt0SdpzRm9pSZLTezqdtsa2BmI8zsq8TXPaPMWBIza2lmI81suplNM7MbEtuzJX9NM/vCzCYl8t+R2L6vmY1JfIYGmVn1qLOWxMxyzOxLMxuWuJ9N2ReY2RQzm2hm4xLbsuKzA2Bm9c1ssJnNNLMZZnZUpudXURD+0gAPA6cD7YFeZtY+2lSlegbovsO2W4D33b0t8H7ifiYqAH7j7u2BI4FfJv57Z0v+zcBJ7n4Y0BHobmZHAvcA97t7G2AVcFV0EUt1AzBjm/vZlB2gm7t33Ob6g2z57AA8CLzj7u2Awwj/HzI7v7tX+htwFDB8m/u3ArdGnasMuVsDU7e5PwvYO/HnvYFZUWcs48/xBnBqNuYHagETgK6Eq2urFveZyqQb0ILwy+gkYBhg2ZI9kW8B0GiHbVnx2QFygfkkTiTKlvzaowiaA4u2ub84sS3bNHH3JYk/LwWaRBmmLMysNdAJGEMW5U8cupkILAdGAHOB1e5ekPiWTP4MPQD8L1CUuN+Q7MkO4MC7ZjbezPomtmXLZ2df4Dvg6cShvyfMrDYZnl9FEVMe/mmS0ec+m1kd4BXgRnf/YdvHMj2/uxe6e0fCv86PANpFm6hszOwsYLm7j486y2441t07Ew4V/9LMjt/2wQz/7FQFOgOPunsnYD07HGbKxPwqiuAboOU291sktmWbZWa2N0Di6/KI85TIzKoRSuI/7v5qYnPW5N/K3VcDIwmHa+qbWdXEQ5n6GToGONvMFgADCYefHiQ7sgPg7t8kvi4HXiMUdbZ8dhYDi919TOL+YEJxZHR+FUUwFmibOPOjOnAxMCTiTOUxBLg88efLCcf+M46ZGfAkMMPd79vmoWzJ39jM6if+vAdhfGUGoTDOT3xbRuZ391vdvYW7tyZ8zj9w90vJguwAZlbbzOpu/TNwGjCVLPnsuPtSYJGZHZjYdDIwnUzPH/UgSabcgDOA2YRjzb+LOk8Z8g4AlgD5hH+lXEU41vw+8BXwHtAg6pwlZD+WsGs9GZiYuJ2RRfk7AF8m8k8F/pjYvh/wBTAHeBmoEXXWUn6OE4Fh2ZQ9kXNS4jZt69/VbPnsJLJ2BMYlPj+vA3tmen5N4SEiIknp0JOIiCSlohARkaRUFCIikpSKQkREklJRiIhIUioKkd1gZi3M7I3ErJ/zzKyfmdUo5TnrStj+ZzM7JfHnG82sVjoyi+wqnR4rUk6JCwfHEKZjeDoxC3F/YJ2735DkeevcvU4pr70AyHP3FanMLFIe2qMQKb+TgE3u/jSE+Z+AXwN9zOw6M+u39RvNbJiZnbjN/fsTa1m8b2aNE9ueMbPzzex6oBkw0sxGVuDPI1IsFYVI+R0MbDe5nofJDRcQJn8rSW1gnLsfDIwCbt/hNR4CviWsudAtlYFFykNFIVLxioBBiT+/QJjSRCRjqShEym860GXbDWZWD2gKfM/2f79qJnkdDRRKRlNRiJTf+0AtM+sD/11S959AP8IqZh3NrIqZtSRMhb1VFX6cqfUSYHQxr70WqJuu4CK7QkUhUk4eThk8FzjfzL4i7EUUuftdwCeEspgOPERYLnWr9cARZjaVMCD+52Jevj/wjgazJRPo9FiRFDGzownTv5/r7hNK+36RbKGiEBGRpHToSUREklJRiIhIUioKERFJSkUhIiJJqShERCQpFYWIiCT1/1IhbCkwHtnHAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# time_sw = [time_sw[i] for i in range(70)]\n",
    "plt.plot(list(range(len(time_qiskit))), time_qiskit)\n",
    "plt.plot(list(range(len(time_sw))), time_sw)\n",
    "plt.xlabel(\"Qubit\")\n",
    "plt.ylabel(\"Time (s)\")\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvwUlEQVR4nO3dd3xV9f348dc7i7DCDHsESFguwAiIC3EBglqcOGtVtK1Wv7a1jrbW/qyr1lYrbcWBbVUQFVFkiQOpk6HMhIQwZGUAgQyyk8/vj89FkpDc3CR3nHvyfj4e90HuOeee+0685p3Pen/EGINSSilVn4hQB6CUUsrZNFEopZTyShOFUkoprzRRKKWU8koThVJKKa+iQh1AIHTt2tUkJCSEOgyllAora9euPWCMia993JWJIiEhgTVr1oQ6DKWUCisi8n1dx7XrSSmllFeuShQiMlVEZuXl5YU6FKWUcg1XJQpjzEJjzIwOHTqEOhSllHINVyUKbVEopZT/uSpRaItCKaX8z1WJQimllP+5KlFo15NSSvmfqxKFdj0ppVypshxWvwy520Py9q5KFEop5Upf/h0W3QvPj4Yl98ORgzXPH94FK56AmWOg+LDf396VK7OVUso1Du2Ez56CpAuhfU9Y9QKsex3O/D/o0BfWvQbbP7PXDhwPRQehdUe/huCqRCEiU4GpiYmJoQ5FKaWazxhY9CuIiIQpf4MOvWHsz+CjP8DHj9hrOvaD8Q/AiOn26wAQN26FmpycbLTWk1Iq7G1+F976MVz0OJz+s5rn9qyFimLoNw4i/DOKICJrjTHJtY+7qkWhlFKuUZJnxyN6nAyjZxx/vs+pQQtFE4VSSoVSVRUs/x1Et4aki6D3KNvV9MmfoDAbpr8BkaH9Va2JQinlLsZA6kJIOBPadA51NLBuDrTpAoMvrPv8xrfgq+ft1yv/DG26woCzIWUBnHYr9A5ey6E+rpoeqwvulFLs+hrm3QBvXg+VFaGNJXc7vPdzeOsmyN1x/PmyIjso3Wsk3LcDLn8ZBp0L2z6B9r3gvN8FP+Y6uCpR6II7pRSb54NEwPdfwKePhjaWlU9DRBRIJCy827Z2qvtqJuTvhYses62fk66Ay1+CX2+Du9ZCrDN+l7kqUbQ0RWUVVFW5b9aaUk1WVQkp78HQKTDqJvj8r5C2NDSxHNwG6+dC8k/gwj/Cjs/gu9eOnS/IsvENuwT6j6v52sgoiI4NbrxeaKIIY08u2cLpT3yMG6c4K9Uk339hB4BPnAaTnoIeJ8G7t8OhOnf4DKyVf4bIaDjzHhj1Y+h/Jix7CPIz7flPHoXKMrjgkeDH1kiaKMJYamYBfTq1QURCHYpSzrBpPkS3tbOHomPhqv+AqbJrESrKghfHgQzY8KYdjG7fw65zuOQ5qCyFRb+ErI22dTHmdug8MHhxNZHOegpTxhhSs/K5dESvUIeilDNUlkPq+zBkIsS0scc6D4RLZ9rB7ZfOg1ZxUF4E5cX2r/2Jj9vZUd5UVcLub2DLIti+Asb+FEZe7/01K5+CyFZwxt3HjnUZBOc+CMt/D5nroXUnOPvXzfqWg0VbFGFq7+FiCkoqGNojLtShKOUMOz6zdY5OmFbz+PBL4MJH7aAy2KmqXZOgtAD+c5mdvlqXzPV2xtLTSTB7EqyaZV/z3p2wYV79cRzYaqe8jr4V2nWreW7sz6HnCMjfY8tu+LkmU6C4qkXRkmo9bcksAGBYT00USgGw6V3bYkg8//hz4+6yj+qKD8G8G2HBHXAwA859yHYRHfrejh9snGfvN/giGHqxvW9EFLx+Jbx7h10gN2zq8e/12ZMQFQvj7j7+XGQUXDkbNr4NyTf75/sOAlclCmPMQmBhcnLybaGOJdBSM/MBGNKjfYgjUcoBKspgy0IYMtn32UKtO8H18+GD/4P/PQ252yCut205SASc9UvbdVR7iur0OfDfH8FbN8P0uZB0vp32uvdbWP+GTQJn/ALaxdf9vp0Hwjn3Ne/7DTJXJYqWZEtWAf27tKFdK/1PqMKUMbB3ra1lFBXTvHtt+8TWRjpxWsPXVhcZDZf83XZFLX8YRGDEtTD+QVuptS6t2sN1b8G/p8Kb18GYOyBtCRxIsy2Jk6+yScZF9LdMmErNzGeotiZUOFs/13b7nPAjuPyV5lVA3TwfYjvCwHMb/1oR23LoO9a2HroNbfg1rTvBDQtg9mT44m/QdwxMfdZ+Lw5ZJOdPmijCUHFZJTsOHuESnfGkwlVBNiy9H9p2s6W043rDRX9q2r3KS2DLYjjhsua1TPqNadz1bbvCrR/ZsY5O/Zv+vmFAZz2FobTsAoxBZzyp8LX4V3aK6s2LbQntr56Hr//VtHtlLIeygsZ3O/lDbJzrkwRoiyIsbfEMZA/XGU8qHKW8Z9c7nPewHRuY+ATk77MtjLhedjprY+/XujMknB2YeJW2KMJRamY+bWMi6dOpdahDUapxinLt1p49T4Fxv7DHIiJh2ovQJxnm32Z/8Rfl+na/ijJI/9DOdgrxng1upj/ZMJSaVcDQnnFERGjpDhVmlj0Ixblww/yav9hj2sD0N+HlC+zaBrBltrufAAlnwOl31Z0Idq6E0jwYNiU48bdQjk8UIjIQeAjoYIy5ItTxhJoxhtTMfC45RQeyVZjJ+AjWz7FlK3qcdPz5tl3g9pW2XEb2ZvvIXA8f/QE6JdgZRbWlfmBrOzVltpPyWUi6nkTkFRHJEZFNtY5PFJE0EckQkfsBjDHbjTG3hCJOJ9qXV0JBSYWuyFbh58u/Q8d+3usbtWoHiefZBWvTXoCffgEd+sGa2cdfW1UFaYvtgjcHleR2o1CNUbwKTKx+QEQigZnAJGA4MF1Ehgc/NGdL3WcHsof11DUUKowU5sCOlXDy1RDVyvfXRUTCqTfaOk4Ht9U8t2e1LSk+tI4yGsqvQpIojDErgdqjVaOBDE8LogyYC1zq6z1FZIaIrBGRNfv37/djtM6yJeto6Q5tUagwsnmBLfd94uWNf+3IG2yNpbWv1jy+ZSFERNe/F7XyGyfNeuoN7K72fA/QW0S6iMi/gJEi8kB9LzbGzDLGJBtjkuPj66mx4gKpmQX066ylO1SY2fQOdBsO3YY1/rXte8CQSbDudagotceMseMTA8525Upop3FSoqiTMeagMeYOY8wgY8zj3q4VkakiMisvLy9Y4QVdala+djup8HJ4N+z+ummtiaNOvdmWEE9daJ/npMChHTrbKUiclCj2An2rPe/jOeYzY8xCY8yMDh3c+RdGcVklOw8c0RXZKrxsnm//bc7K6YHn2plPRwe1tywCBIZc3NzolA+clChWA0kiMkBEYoBrgPcbcwO3tyjSswuoMroHhQozm96B3qc2b8vPiAgYdRN8/znsT7cti76joX13/8Wp6hWq6bFzgK+AISKyR0RuMcZUAHcCy4BUYJ4xZnNj7uv2FsXRPSi060mFjQMZdi1Ec7qdjhp5vR3U/uSPkLUBhmq3U7CEZETUGDO9nuOLgcVBDidsbMkqoG1MJH07tQl1KEr5ZtM7gNS9WK6x2nWzySFlgX2u4xNB46Sup2Zze9dTSmY+Q3q019IdKjwYA5vehv5n2GJ//nB0+9BuJzSvK0s1iqsShZu7nowxbMnM1/EJ5VylhTWfZ2+CA+n+Lf+dcDYkXQhjZvjvnqpBrpqMLyJTgamJiYmhDsXvsvNLyS+p0D2ylTOtfgkW/dJuazp4on2kvgcSCcN9XjfbsIgIuw2pCiptUYSJjBz711pifLsQR6JULcWH4JNHofuJENMW/vc0vDQBvngWBp1rd4JTYc1VLQo3y8gpACCxmyYKFWA7v7AL2kbd6FtdppVPQ/FhuGmhrQp75KCtFLtjpb2HCnuuShRu7nrK2F9I+9go4ts3oqCaUo2VuhDeuhmqyuHrf9rd57zVUsrdDt+8ACOvO1Y6vG0XOOVq+1CuoF1PYSIjp5DEbu0Q0RlPKkA2zIN5N0GvEXDVf0AE3rgSXr/q+MqtR330B4iMhnN/G8xIVZC5KlG4WUbOER2fUIGz9lWYPwP6j4MbFtgB6J9+BRf8P/j+S5g5xo5DlJcce82ur+22pWfcA3E9QxS4CgZNFGEgr6icA4WlOj6hvKsoha0fQVWl768xBr6aCQvvhsTz7YyiVp7PWVSM3UDorrV2iuvKP8M/T4dtn9pNg5Y9CO17wrg7A/P9KMdwVaJw64K7jP06kK0aUFkBb/8EXr/c7iTniyMHYd4N9hf+sKlwzesQ3fr469p3h2mzbEsD4L+XwexJsHctTPidnemkXM1VicKtYxQ/TI3VRKHqUlUF798FWz6wq5VXPF7/mMJRaUvhH2MhfRmc/whc+e+GZzgNOtd2R519n00SPU6GU67x3/ehHMtVicKtMnIKiYmKoI/WeFK1GWNbBOvfgPEPwo8XQ2QMfHCPPVdbaaFNKnOutrWTbvsUzrzHbjnqi+hYmPAQ3L0ebnzP99epsKaJIgxk5BQysGtbIrXGk6rtsyfhm3/C2J/BOffZQeULHrFrGNa9XvPa/H0weyJ89xqceS/c9gn0OLFp79uhN7Tp3Pz4VVjQRBEGMvYXareTOt6aV2w304jr4MI/2emsAKN+DP3GwbKHoDDHHsvaCC+eB7k74dq34PyHfVtMpxQuSxRuHMwuKa9kz6FikrppjSdVTWEOfPg7GDQBpj5nayAdFREBU5+F8iJY8hvYuhxemWgTyU+WQtL5oYtbhSVXJQo3DmZv21+IMTqQrWr59E9QUQKTn4bIOgosxA+2g86b58MbV0HnAXDrR03valItmqtKeLiRzngKQ0W5tqtn/xZbM2l/GrTpAqfdCgPHH+siaqrsFPj2PzD6dugyqP7rzrjb1lxq2xV+9C9opa1S1TSaKBxuW04hEQIJXXXGU1jI3QEvnA2ldttaWneC+KF2FfOWD+zXo2fYaaVNXX/w4W+hVZwdvPYmKgZuWda091CqGk0UDpexv5D+XdrSKkqnIYaFDW9CaQFMfxN6jbRTUEVs6YvN79oZSovuhY8esaudR1wLfU47vpVRVQlVFccPOG/9CLZ9DBc9prOOVNBoonC4jJxCBmmNp/BgDGx8CxLOhCETa56LjoUR021LYvc3sPplWD8X1s6GLomeFkY7uytcdgrkpNrXjZkB4+62FVkrK+DDh+yiutNuC/73p1osTRQOVlFZxY4DR5gwtHuoQ1G+yFwPBzNg3F31XyMC/cbaR0m+Laq3fo4tuAd2LKP7CXZv6MIc+OI5m1TG/szWYNq/Ba5+zXYrKRUkrkoUbtuPYlduEeWVRgeyw8XGtyAiGoZd4tv1sXEw6gb7yNsLEVHHuqqOOvtXdq3Eyqfs837jYOgU/8eulBc6PdbBdMZTGKmqhE3vQNIFTRs76NDbFt+rPVbRbZjdG+L2lXDqzXZ9hO5JooLMVS0Kt8nYbxPFoHitzul4338JBZlw4qOBuX/PU2Dq3wJzb6Ua4KoWhdtkZBfSIy6W9rHRoQ5FNWTjWxDdFoZMCnUkSvmdJgoH0xpPYaKizA5KD71Y92ZQrqSJwqGMMWzL0UQRFrZ9DCWH4aQrQx2JUgGhicKhMvNKOFJWySBNFM638S1o3dlu7KOUC2micKgfZjzpYrvQqmvzn+pKC2HLYjjhMojUsSTlTpooHEqnxoZYVaUt4/2XIbaoX322LIKKYu12Uq7m+EQhIm1F5N8i8qKIXBfqeIJla04hHdtE07WdrsANuqJceP0K+PI5KD4M79xqB6xrO3IAPvqDLcHRd2ywo1QqaEKSKETkFRHJEZFNtY5PFJE0EckQkfs9h6cBbxtjbgN8XPIa/rZmFzC4W3tEF1cFV3YKvHgu7Pzcbgh05WzI2gArHqt5XVUlzL8Nig7CFa/U3DhIKZcJ1af7VaBG1TQRiQRmApOA4cB0ERkO9AF2ey6rDGKMIWOMYWtOIYndtdspqNKWwkvn20qvP14Ep95kp7yOugk+/xvs/OLYtSv/DNs+gclP2cVwSrlYSBKFMWYlkFvr8Gggwxiz3RhTBswFLgX2YJMFeIlXRGaIyBoRWbN///5AhB00+wtLySsuZ7COTwRPVRV8cI/dCW7GCug7+ti5ix6zx9+9HUrybIJY8QScfI1NIkq5nJPay7051nIAmyB6A/OBy0Xkn8DC+l5sjJlljEk2xiTHx8cHNtIA25ptB7KTuuuOZEGTuc6W4Dj9TojrWfNcq3Yw7UXI3wfv3gHv3GY3IJryjNZdUi2C42s9GWOOADf7cq1bqsduzS4AIEm7noInbQlIBCRdWPf5Pslwzm/sWEVMO1uoT1dhqxbCSYliL9C32vM+nmM+M8YsBBYmJyeH9a4u6TmFdGgdTXy7Vg1frPwjbYmdudS2S/3XnPVLKDpgk0n84ODFplSIOanraTWQJCIDRCQGuAZ4vzE3EJGpIjIrLy8vIAEGS0Z2IUnd2umMp2A5vAuyNzZc0C8yCib/2ZYSV6oFCdX02DnAV8AQEdkjIrcYYyqAO4FlQCowzxizuTH3dcN+FMYY0nMKdHwimNKW2n+HTA5tHEo5VEi6nowx0+s5vhhY3NT7umGM4kBhGYeLyknSGU/Bk7YYuiRB1/D93CgVSE7qemo2N7Qojg5kD9YWhX8ZYxfJ1VaSZxfX6T4SStXLVYnCDbbmHJ0aqy0Kv6mqgrd+DDPH2JIc1WV8DFXl2u2klBeuShRuGMxOzy4gLjaKbu11xpPffPL/IGUBHMywi+qqV4RNW2JLhFdfYKeUqsFVicIVXU85hSR11xpPfrNhHnz+DJx6M0z4LWx+F757zZ6rLIety2DwRIiIDG2cSjmYk9ZRtHjGGLZmFzDxxB6hDsUd9qyB9+6EhLPstFaJgB2fwZL7oN9YKMiyYxQ6PqGUV65qUYR719PBI2UcKionsZsOZDdb3l6Yey207wFX/ttuKhQRCT+aBVGx8PbNtjsqMgYGTQh1tEo5mqsSRbh3PaX/MONJB7KbpbzEJomyI3DtmzVXW8f1hMv+AVkbYfVLMOAcW8tJKVUvVyWKcHd0V7skbVE0z7IHbJG/aS9Ct2HHnx8yCUbf7vl64vHnlVI1NDhGISKxwBTgLKAXUAxsAhY1duW08i49u4D2raLoHqcznpps03xY8wqMuwuGepnyesEfbb2mU64NXmxKhSmviUJEHsEmiRXAN0AOEAsMBp7wJJFfGmM2BDhOn4T7yuyt2YUkddcaT012cBu8/wvocxqc97D3a6Nj4bRbgxOXUmGuoRbFKmNMff/HPSMi3YB+fo6pycK9euzWnEIuGNY91GEEzoZ58N1/4dp5EN3av/euKLUD1BGRdmvSyGj/3l+pFszrGIUxZlHtYyISISJxnvM5xpg1gQquJTlYWErukTL3rsg+chAW/wp2rIQ1s/1//w9/C5nr4bJ/QkfH/O2ilCv4NJgtIm+ISJyItMWOT6SIyK8DG1rLku72Xe1WPAalhdBtOHz+Vygr8s99jYGvZsKqWTD2597HJZRSTeLrrKfhxph84DJgCTAAuCFQQbVEGTkunhqbvdkOMJ92C1z8FziSA2t9bFVs/wz+fip89pSd7lpd8WGYdwMse9DWajr/D/6OXCmF74kiWkSisYnifWNMOWC8vyT4wnnBXXp2Ie1bRdEjLjbUofiXMbDkNxDbAcY/AP3H2bULn/+t4VZF1kaYex0UHYRP/wTPjYTVL9vSG3vXwgtn21pNFz4K17wBUTFB+ZaUaml8TRQvADuBtsBKEekP5AcqqKYK5wV3W3MKSHTjjKctH8DO/8G5D0GbzvbY+Adsq2LNy/W/7vAueO0KiI2DO76AW5ZD54Gw6F54PhlevghMFdy81E6FddvPTSkH8SlRGGOeM8b0NsZMNsYYYBdwbmBDazmMMaRlFTDYbQvtyktg2UMQP8wW5Tuq/+kwcLynVXHk+NcV5cJrl0N5MVz3NnTobau73rwEps+FmPZ2odztK6HvacH6bpRqsbwmChG5XkSOu8ZYFSIySETODFx4LUNWfgmHisoZ3isu1KH419cz4fD3MOkJu990deMfhKIDtoxGdeXFMGc6HNoJ09+A7sOPnROxq6p/+jlc/dqxFopSKqAaWkfRBfhORNYCa4H92AV3icA5wAHg/oBG2AKk7LO9eK5KFKUFtsUwZLJtPdTWb4wtxvfFs9Chrx3wzt4E+9ZBYTZcORsS9G8QpZzAa6IwxjwrIs8DE4AzgJOxJTxSgRuMMbsCH6L7pWbaRDG0h4u6nja8CaX5cOa99V8z/gF4+QK7UE4iIX4IDDgbhl8Kw6YEL1allFcN1noyxlQCyz0PFQApmfn079KG9rEuWU1sDKx6EXqOgD7J9V/Xd7QdjI6OteMY0S6b8aWUS7hq46JwrfWUsi+f4T1d1O2083+wfwtcOrPh2Uj9Tw9OTEqpJnNVmfFwnB5bWFrBzoNFDHNTolj1IrTuBCdeHupIlFJ+4KpEEY7SsjwD2U5PFOXFvl2Xtxe2LIJRN/q/8J9SKiR8rfXUXUReFpElnufDReSWwIbWMoTFjKfc7fBkgk0ADVk72y6ES/5JwMNSSgWHry2KV4Fl2I2LANKBewIQT4uTkplPxzbR9Ozg4IHcdXOgosSWCfemohTWvgqDJ0KnhGBEppQKAl8TRVdjzDygCsAYUwFUBiyqFuToQLZjS3dUVcH6ufbrjI9sMqhPyntwZD+MDsvtQJRS9fA1URwRkS54CgGKyFgg/CrvOUxFZRVbsgqcPZD9/ReQtwtOugrKCu1+EvVZ9SJ0HgQDtbqLUm7ia6K4F3gfGCQiXwD/Ae4KWFQtxM6DRyitqHL2QPZ6T22lyU9BTLv6xyn2rYM9q2xrIkLnSCjlJr4WBfwWW7JjHHA7cIJT9skOZ5udPpBddgRSFsAJl9rpronnQdpi2x1V29f/sInklOlBD1MpFVi+znqKBCYD5wEXAneJiJfaDP4jIgM9M67eDsb7BVNKZj4xkREMinfoZkVbFtnupqO//IdOsXWY9n1b87r8fbDpHRh5A7TuGPQwlVKB5WsfwULgx9gige2rPbwSkVdEJEdENtU6PlFE0kQkQ0S8FhU0xmw3xrhyKm7KvnySurcjJsqhXTXr3rD7T/cbZ58nXWBrMtXuflr1op0SO+b24MeolAo4X0t49DHGnNyE+78KPI8d0wB+aJ3MBC4A9gCrReR9IBJ4vNbrf2KMyWnC+4aF1MwCxg+JD3UYdcvbC9tXwDn3HRtzaN0JEs6wieL8h+2xsiN2m9OhU6DzgJCFq5QKHF//lF0iIhc29ubGmJVAbq3Do4EMT0uhDJgLXGqM2WiMmVLr4XOSEJEZIrJGRNbs37+/saEGXU5BCQcKS507kL1xHmDg5KtrHh86BQ6kwcFt9vm6N6DkMJx+Z7AjVEoFia+J4mvgXREpFpF8ESkQkaZuhdob2F3t+R7PsTqJSBcR+RcwUkQeqO86Y8wsY0yyMSY5Pt6hf6VX4+gV2cbYRXZ9x0KXQTXPDZlk/92yyA5qf/0P6J1sK8EqpVzJ166nZ4DTgY2erVCDxhhzELjDl2vDqXpsimcPCkeuodj3rW01TH32+HMd+0GPk22i6JJoy3tc8Tvds1opF/O1RbEb2OSnJLEX6FvteR/PsWYLp+qxKfvy6dOpNR1aO2gPivIS+OYFmHsdRLWG4ZfVfd3Qi2H3N7DiMbs73bBLghqmUiq4fG1RbAdWeIoC/lDDwRjzTBPeczWQJCIDsAniGuDaJtznOOHUokjNdNAeFOXFtkbT53+Dwizofwac93D9U12HTIYVj0PWRrjw0eP3w1ZKuYqvLYodwMdADI2bHjsH+AoYIiJ7ROQWT52oO7FFBlOBecaYzU0JvrZwaVEUlVWw/cARZ3Q77V4Nz58GS++3XUk3fQA3L7Z7Wtenx0nQoZ9dYDfqxuDFqpQKCZ/+FDTGPNKUmxtj6lyma4xZDCxuyj29CZcWRVpWAcaEeHzi6Halyx6EuJ5w00K7X7UvRODiv9iKsrHOTspKqebzmihE5HljzJ0ishBPQcDqjDGO6pw2xiwEFiYnJzu6fGl6dgEAw3o22CgLjNJCWPgLu5o66SKY9oJdI9EYgxs9W1opFaYaalHciO0mejoIsbQY6dmFxEZH0LdTm+C/+YGtdrD64FaY8Ds4814t4qeU8qqhRLENwBjzWRBiabZw6XpKzy4gsVs7IiKCPKV0y2KYPwOiYuD6+TBIy4ErpRrWUKKI91b8r4mzngImXLqetmYXMm5Ql+C9YVUVfPYkfPYE9BwBV78GHfs2+DKllIKGE0Uk0A7Q1VR+kldcTlZ+CUndAzQ+kbURKssguo19iMDi+yB9CZxyLUx5BqJbB+a9lVKu1FCiyDTG/DEokfhBOHQ9ZeTYgezB3QNQWvz7L2H2pOOPR0TB5KfhtFt1BbVSqtEaShRh9VslHLqe0rMLARgciBbF2n9DqziYNssuoisvgrIi6Hsa9Brp//dTSrUIDSWK84ISRQuSnl1A6+hIenf0c/dPSR6kvAcjph8r3KeUUn7gdV6kMaZ2iXDVTOnZBSR1D8CMp03zoaIYRl7v3/sqpVo8V02gF5GpIjIrLy8v1KHUKz27kKRuAeh2+u6/0G049Brl/3srpVo0VyUKp9d6OlxUxv6CUv8PZGenwN61tjWhg9VKKT9zVaJwuh8Gsnv4uUWx7nWIiD5+NzqllPIDTRRBdLTGk19nPFWUwfq5dgC7bVf/3VcppTw0UQTR1uwC2rWKoleHWP/dNH0pFB2AkTf4755KKVWNqxKF0wez07MLSezWDmnKOEJ5Mfw9Gf47DfasPXb8u9egfU8YNMF/gSqlVDWuShROH8zemlPQ9IHsjI9txdddX8FLE+CNqyH9Q8hYDiOu1V3mlFIBo79dgiT3SBkHCsuaPj6RssDuGfGL72D1S/Dl3223E8CI6/wWp1JK1aaJIkiODmQ3qRhgeQmkLYUTLrPJ4uxfw+gZ8PW/wFRBl0H+DVYpparRRBEkW7ObUQxw28dQVmATxVGxHWD8b/wTnFJKeeGqMQonS88upH2rKHrENWHG0+YFtiUx4By/x6WUUg3RRBEkR2s8NXrGU0UppC2BoRdDZHRgglNKKS9clSicOj3WGEN6dkHTBrK3fWK7nYb/yP+BKaWUD1yVKJw6PfZAYRmHisqbNpC9eQHEdoSB2u2klAoNVyUKp2ryQHZFKaQthqFTtNtJKRUymiiCoMk1nrZ9CqX5NWc7KaVUkGmiCIL0nELiYqPo1r5V416YssBOg9XZTkqpENJEEQRbMvMZ2jOucTOeKkphi6fbKSomcMEppVQDdMFdgFVWGVIzC7hmdF/vF5bkwb51kLcbDu+C7M1QmgfDLwtGmEopVS9NFAG248ARissrGd4zrv6LjhyEF86G/D2eAwJxvWDwJBg4PhhhKqVUvTRRBNjmfXZNxwm96pmyawy893M4kgNX/Rd6nARxvbW7SSnlGI5PFCJyGXAxEAe8bIz5MLQRNU5KZj4xkREkdqtnauw3L0D6Epj4JAy/JLjBKaWUDwI6mC0ir4hIjohsqnV8ooikiUiGiNzv7R7GmAXGmNuAO4Cw2xQ6ZV8+Sd3bERNVx486cz0s/53tYhpze/CDU0opHwR61tOrwMTqB0QkEpgJTAKGA9NFZLiInCQiH9R6dKv20t96Xhc2jDFs3pfPCb3qGJ8oLYS3fwJtusClM6Epu94ppVQQBLTryRizUkQSah0eDWQYY7YDiMhc4FJjzOPAlNr3EDun9AlgiTHm2/reS0RmADMA+vXr559voJmy8kvIPVJW9/jEkvvg4Da4aSG07RL84JRSykehWEfRG9hd7fkez7H63AWcD1whInfUd5ExZpYxJtkYkxwfH++fSJspZV8+wPEtik3zYd3rcPavYMBZIYhMKaV85/jBbGPMc8BzvlwrIlOBqYmJiYENykeb9+UjAkOrT40tyIJF90LvU+Ecr8MzSinlCKFoUewFqq8+6+M51mxOqx67eV8eCV3a0q6VJx8bAwvvhvJiuOxfEOn4PK2UUiFJFKuBJBEZICIxwDXA+/64sdP2o9i8L5/h1budvnsN0pfCeQ9D/ODQBaaUUo0Q6Omxc4CvgCEiskdEbjHGVAB3AsuAVGCeMWazP97PSS2KvOJy9hwqPrYi+9D3sPQBSDgLxtQ71KKUUo4T6FlP0+s5vhhY7O/3c9IYRY2B7Koqu/oaY6fCRmgtRqVU+HDVbywntShqlO74+h+w839w0WPQqX+II1NKqcZxVaIIio8ega9m2oFpL1L25dOtfSvi9yy3q6+HToFRNwYpSKWU8h9XTbsJeNdTYQ58/oz9OnMDXPIcRNW9GVFKZj6XdN4N79wPvUbBtFm6+lopFZZc1aIIeNdTxkf231Omw4a58J9LbYnwWkrKK6nMSeOXB39vK8Fe+ybEtA1MTEopFWCuShQBt3U5tOsOl/4DrngF9n0HL02A7JQaXVHbt2cwO/oJIiJj4Pp3oG3XEAatlFLNo11PvqqsgG0fw9CpdtbSiZdDx/4wZzr883SIjIHWnaFNZxLyD1BFIQcvXUCvzgP8H4tSSgWRqxKFMWYhsDA5Ofk2v998z2q7XWnSBceO9UmGGStg0ztQdACKcqH4EHuL2/F48UReGjrW72EopVSwuSpRBNTWD0Eij9+atENvOOMXNQ795h9fENU2gogIHbxWSoU/HaPwVcZy6DcWWnf0ellllSE1s6Bm6Q6llApjrkoUfqn1VFV5/LH8fZC1sWa3Uz3SswsoLq/kpN6hX/SnlFL+4KpE0azpscbA4vs8pTZqOTotNrHhRLF6Zy4Aowd0bnwMSinlQK5KFM0iAq07wfo5sHlBzXNbl0P7XtD9hAZvs2pHLj3iYunTqXVg4lRKqSDTRFHd2b+yq6g/uAfyM+2xynLY9qntdmpgZbUxhtU7cxk9oDOiq7CVUi6hiaK6yGiY9iKUl8B7P7PdUbu+hrICn8YnduUWkZ1fymna7aSUchFXJQq/DGZ3TYSLHoVtn8CqF+202IhoGHBOgy9dtcMzPpGgiUIp5R6uShR+q/WUfAskXWirvm56B/qfDrENT3ddvTOXDq2jSerWrnnvr5RSDuKqROE3InDJ8xDdBvL32qThg9U7D3FaQmddaKeUchVNFPVp393uRte6Ewy9uMHLcwpK2HHgCKMHdApCcEopFTxawsOboZNhyA6f9pFYveMQAKfp+IRSymW0RdEQH6e5rt6ZS+voSE7UFdlKKZfRROEnq3bkMqp/R6Ij9UeqlHIXV/1W88v02CbILyknNStfu52UUq7kqkQR8K1Q67F25yGM0fUTSil3clWiCJVVO3OJihBG9tMZT0op99FE4Qerd+RyUp8OtI6JDHUoSinld5oomqmkvJINe/K020kp5VqaKJpp3e7DlFVW6f4TSinX0kTRTCvT9yMCyf01USil3EkTRTOUVlQyb81uzhvajQ5tokMdjlJKBYQmimZYtCGTA4Vl3DQuIdShKKVUwDg+UYjIMBH5l4i8LSI/DXU8RxljePXLnQyKb8uZiV1DHY5SSgVMQBOFiLwiIjkisqnW8YkikiYiGSJyv7d7GGNSjTF3AFcBZwQy3sb4bvdhNuzJ46ZxCbrtqVLK1QLdongVmFj9gIhEAjOBScBwYLqIDBeRk0Tkg1qPbp7XXAIsAhYHOF6f/fvLnbRvFcW0UX1CHYpSSgVUQMuMG2NWikhCrcOjgQxjzHYAEZkLXGqMeRyYUs993gfeF5FFwBt1XSMiM4AZAP369fPPN1CPnIISFm/M5Pqx/WnXSiu1K6XcLRS/5XoDu6s93wOMqe9iERkPTANa4aVFYYyZBcwCSE5ONn6Is15vfLOL8krDjacnBPJtlFLKERz/57AxZgWwwpdrRWQqMDUxMTFg8ZRVVPH6N7sYPySeAV3bBux9lFLKKUIx62kv0Lfa8z6eY80WjOqxSzZlsr+gVKfEKqVajFAkitVAkogMEJEY4BrgfX/cOND7UVRVGV75fAcDurblnKT4gLyHUko5TaCnx84BvgKGiMgeEbnFGFMB3AksA1KBecaYzf54v0C3KOau3s36PXn8dPwgIiJ0SqxSqmUI9Kyn6fUcX0wAproGcowiO7+Ex5ekcvrALlx5qk6JVUq1HI5fmd0YgWxRPPzeZsoqqnhs2km6wE4p1aK4KlEEytJNWSzdnMU95w/WmU5KqRbHVYkiEIPZecXl/P69TQzrGcetZw3w232VUipcuCpRBKLr6cmlWzhQWMqTl59EdKSrflxKKeUTxy+4C6YH5m/gmx25VFQaKiqrqKgy5BSUcttZAzi5T8dQh6eUUiHhqkTR3FlPfTq1oaCkgqgIISoyguhIIb5dK346PnArvZVSyunEmICWRQqJ5ORks2bNmlCHoZRSYUVE1hpjkmsf1053pZRSXmmiUEop5ZWrEkWgaz0ppVRL5KpEEYzqsUop1dK4KlEopZTyP00USimlvNJEoZRSyitXJQodzFZKKf9z5YI7EdkPfN/El3cFDvgxnEDTeANL4w28cIvZzfH2N8Yct32nKxNFc4jImrpWJjqVxhtYGm/ghVvMLTFeV3U9KaWU8j9NFEoppbzSRHG8WaEOoJE03sDSeAMv3GJucfHqGIVSSimvtEWhlFLKK00USimlvNJEUY2ITBSRNBHJEJH7Qx1PbSLyiojkiMimasc6i8hyEdnq+bdTKGOsTkT6isinIpIiIptF5G7PcUfGLCKxIrJKRNZ74n3Ec3yAiHzj+Vy8KSIxoY61OhGJFJHvROQDz3PHxisiO0Vko4isE5E1nmOO/DwAiEhHEXlbRLaISKqInO7UeEVkiOfnevSRLyL3+CNeTRQeIhIJzAQmAcOB6SIyPLRRHedVYGKtY/cDHxtjkoCPPc+dogL4pTFmODAW+LnnZ+rUmEuBCcaYU4ARwEQRGQs8CfzVGJMIHAJuCV2IdbobSK323OnxnmuMGVFtbr9TPw8AzwJLjTFDgVOwP2dHxmuMSfP8XEcApwJFwLv4I15jjD7sgP7pwLJqzx8AHgh1XHXEmQBsqvY8Dejp+bonkBbqGL3E/h5wQTjEDLQBvgXGYFe1RtX1OQn1A+jj+Z9/AvABIA6PdyfQtdYxR34egA7ADjyTfpweb60YLwS+8Fe82qI4pjewu9rzPZ5jTtfdGJPp+ToL6B7KYOojIgnASOAbHByzpxtnHZADLAe2AYeNMRWeS5z2ufgbcB9Q5XneBWfHa4APRWStiMzwHHPq52EAsB+Y7enae0lE2uLceKu7Bpjj+brZ8WqicBFj/2Rw3HxnEWkHvAPcY4zJr37OaTEbYyqNbbr3AUYDQ0MbUf1EZAqQY4xZG+pYGuFMY8wobBfvz0Xk7OonHfZ5iAJGAf80xowEjlCr28Zh8QLgGZO6BHir9rmmxquJ4pi9QN9qz/t4jjldtoj0BPD8mxPieGoQkWhsknjdGDPfc9jRMQMYYw4Dn2K7bjqKSJTnlJM+F2cAl4jITmAutvvpWZwbL8aYvZ5/c7D956Nx7udhD7DHGPON5/nb2MTh1HiPmgR8a4zJ9jxvdryaKI5ZDSR5ZozEYJtu74c4Jl+8D9zk+fom7DiAI4iIAC8DqcaYZ6qdcmTMIhIvIh09X7fGjqekYhPGFZ7LHBOvMeYBY0wfY0wC9vP6iTHmOhwar4i0FZH2R7/G9qNvwqGfB2NMFrBbRIZ4Dp0HpODQeKuZzrFuJ/BHvKEedHHSA5gMpGP7pR8KdTx1xDcHyATKsX/t3ILtk/4Y2Ap8BHQOdZzV4j0T28zdAKzzPCY7NWbgZOA7T7ybgN97jg8EVgEZ2OZ8q1DHWkfs44EPnByvJ671nsfmo/+POfXz4IltBLDG85lYAHRyeLxtgYNAh2rHmh2vlvBQSinllXY9KaWU8koThVJKKa80USillPJKE4VSSimvNFEopZTyShOFUs0gIn1E5D1PZc7tIvK8iLRq4DWF9Rz/o4ic7/n6HhFpE4iYlWosnR6rVBN5FhR+gy3xMNtTgXgWUGiMudvL6wqNMe0auPdOINkYc8CfMSvVFNqiUKrpJgAlxpjZYOtEAf8H3Cgid4rI80cvFJEPRGR8ted/9ex58bGIxHuOvSoiV4jIL4BewKci8mkQvx+l6qSJQqmmOwGoUZDP2KKHO7EF5erTFlhjjDkB+Ax4uNY9ngP2YfdtONefASvVFJoolAq+KuBNz9evYUudKOVYmiiUaroU7E5iPxCROKAHtt5O9f+/Yr3cRwcKlaNpolCq6T4G2ojIjfDDdrp/AZ7H7ow2QkQiRKQvtpz2UREcq+56LfB5HfcuANoHKnClGkMThVJNZOyUwR8BV4jIVmwrosoY8yfgC2yySAGew26retQRYLSIbMIOiP+xjtvPApbqYLZyAp0eq5SfiMg4bCn4Hxljvm3oeqXChSYKpZRSXmnXk1JKKa80USillPJKE4VSSimvNFEopZTyShOFUkoprzRRKKWU8ur/A/tOkLFRQG/1AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "time_sw = [time_sw[i] for i in range(70)]\n",
    "plt.plot(list(range(len(time_qiskit))), time_qiskit)\n",
    "plt.plot(list(range(len(time_sw))), time_sw)\n",
    "plt.xlabel(\"Qubit\")\n",
    "plt.ylabel(\"Time (s)\")\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[4 6 5 7]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([[2, 0],\n",
       "        [3, 0],\n",
       "        [2, 1],\n",
       "        [3, 1]]),\n",
       " array([0, 0, 2, 3]))"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "# np.kron()\n",
    "a = np.array([0, 1])\n",
    "b = np.array([2, 3])\n",
    "# np.tensordot(a, a)\n",
    "\n",
    "def cartesian_coord(arr1, arr2):\n",
    "    grid = np.meshgrid(arr2, arr1)\n",
    "    coord_list = [entry.ravel() for entry in grid]\n",
    "    # print(coord_list[0] << 1 | coord_list[0])\n",
    "    print(coord_list[0] << 1 | coord_list[1])\n",
    "    points = np.vstack(coord_list).T\n",
    "    return points\n",
    "\n",
    "cartesian_coord(a, b), np.kron(a, b)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.82 , 0.14 ],\n",
       "       [0.14 , 0.905]])"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import scipy as sp\n",
    "def nearest_unitary(A):\n",
    "    \"\"\"\n",
    "    Calculate the closest unitary to a given matrix.\n",
    "\n",
    "    Calculate the unitary matrix U that is closest with respect to the\n",
    "    operator norm distance to the general matrix A.\n",
    "\n",
    "    D.M.Reich. \"Characterisation and Identification of Unitary Dynamics\n",
    "    Maps in Terms of Their Action on Density Matrices\"\n",
    "\n",
    "    Args:\n",
    "        A (np.ndarray): The matrix input.\n",
    "\n",
    "    Returns:\n",
    "        (np.ndarray): The unitary matrix closest to A.\n",
    "        Return U as a numpy matrix.\n",
    "\n",
    "    Thank you to Ed Younis, this is based on code from baseline.qfast\n",
    "    \"\"\"\n",
    "    try:\n",
    "        if len(A.shape) == 2 and A.shape[0] != A.shape[1]:\n",
    "            raise TypeError(\"A must be a square matrix.\")\n",
    "\n",
    "        V, __, Wh = sp.linalg.svd(A)\n",
    "        U = np.array(V.dot(Wh))\n",
    "        return U\n",
    "    except Exception:\n",
    "        return A\n",
    "    \n",
    "measure_fids0 = 0.90\n",
    "measure_fids1 = 0.95\n",
    "M = np.array([[measure_fids0, 1 - measure_fids0],\n",
    "    [1 - measure_fids1, measure_fids1]])\n",
    "M @ M.T.conj()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.00000000e+00, 1.64340319e-16],\n",
       "       [1.64340319e-16, 1.00000000e+00]])"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "uM = nearest_unitary(M)\n",
    "uM @ uM.T.conj()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[0.9 , 0.1 ],\n",
       "        [0.05, 0.95]]),\n",
       " array([[ 0.99963497,  0.02701716],\n",
       "        [-0.02701716,  0.99963497]]))"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M, uM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[1, 0, 0, ..., 1, 1, 0],\n",
       "        [1, 0, 0, ..., 1, 1, 0],\n",
       "        [1, 0, 0, ..., 1, 1, 0],\n",
       "        ...,\n",
       "        [0, 0, 1, ..., 0, 1, 0],\n",
       "        [0, 1, 0, ..., 0, 0, 1],\n",
       "        [0, 0, 1, ..., 0, 1, 0]],\n",
       "\n",
       "       [[1, 0, 0, ..., 1, 1, 0],\n",
       "        [1, 0, 0, ..., 1, 1, 0],\n",
       "        [1, 0, 0, ..., 1, 0, 0],\n",
       "        ...,\n",
       "        [0, 1, 0, ..., 0, 1, 1],\n",
       "        [1, 0, 0, ..., 1, 1, 1],\n",
       "        [0, 1, 0, ..., 1, 0, 1]],\n",
       "\n",
       "       [[1, 0, 0, ..., 1, 1, 0],\n",
       "        [1, 0, 0, ..., 0, 1, 0],\n",
       "        [1, 0, 0, ..., 1, 1, 0],\n",
       "        ...,\n",
       "        [1, 0, 1, ..., 0, 1, 1],\n",
       "        [1, 1, 1, ..., 0, 0, 0],\n",
       "        [0, 1, 1, ..., 1, 0, 0]],\n",
       "\n",
       "       ...,\n",
       "\n",
       "       [[1, 0, 0, ..., 1, 1, 0],\n",
       "        [0, 0, 0, ..., 1, 1, 0],\n",
       "        [1, 0, 0, ..., 1, 1, 0],\n",
       "        ...,\n",
       "        [1, 0, 0, ..., 0, 0, 0],\n",
       "        [0, 1, 1, ..., 1, 0, 1],\n",
       "        [1, 0, 1, ..., 1, 0, 1]],\n",
       "\n",
       "       [[1, 0, 0, ..., 1, 1, 0],\n",
       "        [1, 1, 0, ..., 1, 1, 0],\n",
       "        [1, 0, 0, ..., 1, 1, 0],\n",
       "        ...,\n",
       "        [1, 1, 0, ..., 1, 1, 1],\n",
       "        [1, 1, 0, ..., 0, 1, 1],\n",
       "        [1, 1, 1, ..., 0, 0, 0]],\n",
       "\n",
       "       [[1, 0, 0, ..., 1, 1, 0],\n",
       "        [1, 0, 0, ..., 1, 1, 0],\n",
       "        [0, 1, 0, ..., 1, 0, 0],\n",
       "        ...,\n",
       "        [0, 0, 1, ..., 1, 1, 0],\n",
       "        [1, 0, 1, ..., 0, 0, 0],\n",
       "        [1, 0, 0, ..., 1, 1, 1]]], dtype=uint8)"
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "import scipy.io as sio\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import itertools\n",
    "from collections import defaultdict\n",
    "from tqdm import tqdm\n",
    "\n",
    "from jax import numpy as jnp\n",
    "import jax\n",
    "import time\n",
    "from functools import lru_cache\n",
    "from line_profiler import LineProfiler\n",
    "\n",
    "qnum = 24\n",
    "stats01s = sio.loadmat(f'{qnum}qubit_data/stats_0.mat')['stats01s']\n",
    "stats01s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "tf2",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
